This notebook utilizes the data from Martin et al 2019 and processed using Seurat pipeline. Here we explore the data and generate figures to understand the CD ileal data quality and the cell type differences with inflammation. The data used here is availabe to download as rds file through our IBDGC DataCommons

Package installation

Uncomment the installation line if the environment doesnot have those packages installed.

#install.packages(c("ggpubr", "Seurat", "ggplot2", "dplyr","tidyr","viridis"))
#devtools::install_github("rpolicastro/scProportionTest")

#Load the libraries
library(ggpubr)
library(Seurat)
library(ggplot2)
library(dplyr)
library(tidyr)
library(viridis)

Reading data and QC

setwd('/Users/mamtagiri/Downloads/CD ileal seurat')
test<-readRDS("/Users/mamtagiri/Downloads/CD ileal seurat/Immune.rds")
test$annotation<-test@active.ident

##### UMI plot for the clusters ####
meta<-test@meta.data
p1=ggplot(meta, aes(x=nCount_RNA)) + geom_histogram(binwidth=1000, fill="#69b3a2", color="#e9ecef", alpha=0.9) 
p2=ggplot(meta, aes(x=annotation, y=nCount_RNA,fill=seurat_clusters))+geom_boxplot(show.legend = FALSE)+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1,size=15))

ggarrange(p1,p2, labels = c("UMI histogram","UMAP distribution"),common.legend = FALSE,ncol = 1, nrow = 2)

Table with cell number per condition per cluster

We explore how the cell contribution for each cell type is across the two conditions (inflamed, uninflamed)

##### Table with cell number per condition per cluster
table(test@active.ident,test@meta.data$group)
##                           
##                            inflamed uninflamed
##   Epithelial                    102        366
##   Central Trm                  5621       4244
##   CD8 cytotoxic T cells        2665       4441
##   IgA plasma cells             2955       3315
##   gd T cells                   2234       1025
##   Naive B cells                2134       1070
##   CD8+ Trm                      951       1532
##   Trm                          1154       1318
##   Mem. B cells                 2716       1276
##   Type 3 cytokine Trm           598       1223
##   IgM plasma cells              674        932
##   ILC1                          830        651
##   Tregs                        1065        267
##   DC2                           614        687
##   IgG plasma cells             2126        381
##   Residential macrophages       535        582
##   ILC3                          472        634
##   Fibroblasts                   417        731
##   moDCs                         319        445
##   ACKR1+ endothelial            486        219
##   CD36+ endothelial             125        491
##   Inflammatory macrophages      570         23
##   DC1                           335        153
##   Activated fibroblasts         370          3
##   Activated DCs                 315         20
##   Smooth muscle                 132        128
##   Glial cells                    29        220
##   Lymphatics                    107        105
##   Plasmablasts                  133         59
##   Pericytes                      27         96

UMAP of the 11 paired inflamed, uninflamed ileal single cell data

## UMAP of the 11 paired inflamed, uninflamed ileal single cell data ###
p1=DimPlot(test,reduction = "umap",label = TRUE,label.size = 6.5,repel = TRUE)+NoLegend()
p2=DimPlot(test,reduction = "umap",label = TRUE,label.size = 5,repel = TRUE,split.by = "group")+NoLegend()
ggarrange(p1,p2, labels = c("UMAP","UMAP by condition"),common.legend = FALSE,ncol = 1, nrow = 2)

Make a cell proportion bar graph

#### Make a cell proportion bar graph ###
meta<-test@meta.data
counts <- as.data.frame(table(test@active.ident,test@meta.data$group))
ggplot(counts, aes(fill=Var2, y=Freq, x=Var1)) + 
  geom_bar(position="fill", stat="identity",alpha=0.8)+ scale_fill_brewer(palette="Set2") + xlab(" Clusters")+ ylab("Frequency") + labs(fill = "Condition")+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1,size=15))

Perform a permutation test for cell proportion changes

library("scProportionTest")
prop_test <- sc_utils(test)

prop_test <- permutation_test(
  prop_test, cluster_identity = "annotation",
  sample_1 = "uninflamed", sample_2 = "inflamed",
  sample_identity = "group"
)
permutation_plot(prop_test)

Cell Proportion per patient per condition

Visualize how consistent these inflammation related cell type changes are in each patient

counts <- as.data.frame(table(test@active.ident,test@meta.data$group,test$SampleID))

#Renaming cells under broader category
Plasma<-c("IgA plasma cells","IgG plasma cells","IgM plasma cells","Plasmablasts")
MNP<-c("Residential macrophages","Inflammatory macrophages","moDCs","pDCs","DC1","DC2","Activated DCs")
Tc<-c("Central Trm","CD8 cytotoxic T cells","gd T cells","CD8+ Trm","Trm","Type 3 cytokine Trm","Tregs")
Stromal<-c("Smooth muscle","Fibroblasts","Activated fibroblasts","Lymphatics","Pericytes","CD36+ endothelial","ACKR1+ endothelial","Glial cells")
B<-c("Mem. B cells","Naive B cells")
Mast<-c("Mast cells")
ILC<-c("ILC3","ILC1")
Epithelial<-c("Epithelial")

counts$id<-paste0(counts$Var3,"-",counts$Var2)

counts_new<-counts%>%mutate(CellGroup = case_when(
  Var1%in%Plasma ~ "Plasma",
  Var1%in%MNP ~ "MNP",
  Var1%in%Tc ~ "Tcells",
  Var1%in%Stromal ~ "Stromal",
  Var1%in%B ~ "B cells",
  Var1%in%Mast ~ "Mast",
  Var1%in%Epithelial ~ "Epithelial",
  Var1%in%ILC ~ "ILCs"
))

ggplot(counts_new, aes(fill=CellGroup, y=Freq, x=id)) + geom_bar(position="fill", stat="identity",alpha=0.8)+ scale_fill_brewer(palette="Set2") + xlab(" Clusters")+ ylab("Frequency") + labs(fill = "Condition")+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1,size=7))

Overlay Cell type markers from Martin et al. 2019

### Overlay Cell type markers from Martin et al. 2019
featuresimmune<- c("CD3D", "CD2", "CD7", "TNFRSF17", "XBP1", "BANK1", "CD79B", "CD22", "MS4A1", "HLA-DRB1", "HLA-DQA1", "LYZ", "IL3RA", "IRF7", "GZMB", "LILRA4", "CLEC4C", "TPSAB1", "CMA1", "KIT", "FCER1A", "FCER1G", "KLRB1", "TYROBP", "PLVAP", "VWF", "COL3A1", "COL1A1", "ACTA2")
test<-ScaleData(test,features = featuresimmune)
DoHeatmap(subset(test,downsample=400),size=7,angle=90,lines.width=13,features = featuresimmune)+NoLegend() +scale_fill_viridis()+theme(axis.text.y = element_text(size = 15,colour = "black"))